Superionic effect and anisotropic texture in Earth’s inner core driven by geomagnetic field

Seismological observations suggest that Earth’s inner core (IC) is heterogeneous and anisotropic. Increasing seismological observations make the understanding of the mineralogy and mechanism for the complex IC texture extremely challenging, and the driving force for the anisotropic texture remains unclear. Under IC conditions, hydrogen becomes highly diffusive like liquid in the hexagonal-close-packed (hcp) solid Fe lattice, which is known as the superionic state. Here, we reveal that H-ion diffusion in superionic Fe-H alloy is anisotropic with the lowest barrier energy along the c-axis. In the presence of an external electric field, the alignment of the Fe-H lattice with the c-axis pointing to the field direction is energetically favorable. Due to this effect, Fe-H alloys are aligned with the c-axis parallel to the equatorial plane by the diffusion of the north–south dipole geomagnetic field into the inner core. The aligned texture driven by the geomagnetic field presents significant seismic anisotropy, which explains the anisotropic seismic velocities in the IC, suggesting a strong coupling between the IC structure and geomagnetic field.

superionic material Li2O at high temperatures 3 . The elastic properties of crystals are expressed as the relationship between stress and strain: where refers to a stress tensor, refers to a strain tensor, and represents a fourth-order elastic modulus. The elastic constant was calculated by distorting the equilibrium structure and solving the stress-strain relations (see method for details; Table 1).
The statistical errors of the calculated temperature and stress were evaluated using the blocking method 4 . The standard deviation is defined as follows 5 : (5) where n is the number of data points. The standard deviations of the time-averaged properties in the NVT simulations are counted, and the errors are 3-9 K in temperature and 0.02-0.15 GPa in stress. The errors are ~2% in .
It is worth noting that the thermodynamic properties for superionic Fe alloys is not well constrained; thus, the estimated adiabatic elastic constants may present some deviation. In addition, as the thermal conductivity is high in the inner core, the presence of thermal convection, and possible Joule heating due to the penetration of the outer core magnetic field, it is still not clear if the seismic wave velocities in the inner core are isothermal, intermediate or adiabatic. In this case, we still use the isothermal data for discussion. The adiabatic compressional velocity is higher than the isothermal velocity, but the difference in compressional velocity anisotropy is negligible (Supplementary Fig. 1 and Supplementary Table 2).
We calculated the bulk modulus B and shear modulus G using the Voigt average scheme, which is proven to be more appropriate and accurate in calculating the seismic wave properties 9 .
The compressional wave velocity , shear wave velocity , and bulk sound velocity are calculated by: The results are provided in Table S1. We also compared the compressional velocities in superionic FeH alloys with the PREM data 10

Supplementary Discussion 2. Reversal of the fast velocity axis and change in the c/a ratio in superionic Fe-H alloys
In this section, we provide a detailed method for seismic velocity anisotropy calculation and a discussion on the relationship between seismic anisotropy and the lattice c/a ratio.
To further investigate the elastic anisotropy, we calculated the azimuthal angle-dependent velocity of superionic Fe-H alloys by solving the following Christoffel equation 13,14 .
where is the propagation direction and is the polarization direction. Here, we investigated velocity anisotropy by analysing the velocity change as a function of ray angle with Earth's rotation axis. For a hexagonal system with cylindrical symmetry, we have: where is the density, VP is the compressional velocity, and is the angle between the ray path in the IC and Earth's rotation axis. The compressional wave velocity anisotropy and the maximum shear wave splitting anisotropy in Fe alloys (Supplementary Table 1) are calculated by: = [ The velocity anisotropies are also calculated by solving the Christoffel equation employing MTEX 15 , as shown in Supplementary Fig. 3.
Interstitial H impurities have a significant influence on the seismic velocity anisotropy in hcp Fe. The reversal of the fast velocity direction from the c-axis to the basal plane is observed in the Fe-H alloy with increasing temperature and H content (Fig. 1). The c/a ratio is an important parameter governing the elastic anisotropy of the hcp structure [16][17][18][19] . The c/a ratio of hcp-Fe increases with temperature and approaches the ideal value (1.633) at the temperature condition of the IC 20,21 . As a result, the elastic anisotropy of hcp-Fe under IC conditions is supposed to be weak. However, a recent experimental study suggests that the c/a ratio of hcp-Fe is substantially lower than the ideal value, indicating significant anisotropy in hcp-Fe under IC conditions 19 ( Supplementary Fig. 4). We calculated the cell parameters using AIMD simulations at different temperatures and pressures under hydrostatic conditions. Our calculated c/a ratio of pure Fe at 360 GPa and 6000 K is approximately 1.61, which is consistent with previous experimental measurements, and the c/a ratio of the Fe-H alloy is larger than that of pure Fe. The c/a ratio grows further with temperature and H content, but it is still lower than the ideal value. Unlike pure Fe, the seismic anisotropy shows a trend of initially decreasing and then increasing with the c/a ratio. This effect is mainly caused by the significant decrease in C33 with increasing temperature. It is likely that the highly diffusive H-ions taking more unstable interstitial sites at high temperature within the basal plane are responsible for the elastic softening in C33

Supplementary Discussion 3. Elastic anisotropy in Fe-H and other Fe alloys
Here, elastic anisotropies in Fe alloys are compared and discussed. Using the calculated elastic constants of previous computational studies, the compressional wave velocity anisotropies in pure
In considering the velocity change in FeH0.25 with azimuthal angle, we need to deposit its basal plane parallel to Earth's rotation axis, as shown in Supplementary Fig. 7. Here, we considered two models for the configuration of  Supplementary Fig. 7a). In this model, the change in the velocity as a function of angle ξ between the ray path and Earth's rotation axis is calculated by averaging the velocities of all the possible propagation direction paths corresponding to the different alignment of crystal with a rotatable c-axis in the equatorial plane ( Supplementary Fig. 7b). As shown in Supplementary   Fig.7b, the crystal has a rotational degree of freedom in the equatorial plane. The propagation direction vector is ( 1 , 2 , 3 ). Thus, we have: { 1 = cos 2 = sin sin 3 = sin cos (18) and we can get ( , ) by substituting (cos , sin sin , sin cos ) into Eq. 14 and solving the Christoffel equation. In the IEP model, we set the sampling interval of as 1° and ( ) was calculated by: It can be viewed as a cylindrically averaged (averaged over ) aggregate with a basal plane parallel to Earth's rotation axis. These two models are compared with the seismic travel time data, as discussed in the main text (Fig. 4). Here, PKIKP refers to the compressional wave path transmitted through the inner core. PKPab and PKPbc are used as reference phases traversing only the mantle and outer core. Generally, the AEP model presents the slowest angle, which fits the IMIC absolute PKIKP data better 25 . This may be due to insufficient ray path sampling of the IMIC.
Alternatively, the IMIC may present a highly anisotropic equatorial plane, which is also suggested by recent observations of the coda-correlation wavefield. The equatorial anisotropy was proposed based on earthquake coda correlation, which is able to sample the deepest part of the IC 35,36 . As suggested by recent studies 37,38 , earthquake coda correlation can be an efficient approach to provide additional constraints on Earth's interior structure if the real reconstructed ray paths from the coda correlation wavefield are carefully interpreted. Understanding the discrepancy between coda interferometry and ballistic observations 39,40 may further provide critical clues for understanding the IC anisotropy structure. In addition, with the increase in seismological observations and the application of the coda correlation method, more information from the deepest part of Earth will be collected, and the mystery of the IMIC may be revealed. In particular, if the proposed equatorial anisotropy of the IMIC is confirmed with additional data, it should be critical evidence of the presence of a superionic Fe-H alloy in the IC.
We also compared the velocity anisotropy of other Fe alloys (360 GPa) 5,11,12 with the IMIC observation data 25  However, the ray paths for these waves are different in the mantle; thus, the residuals are also affected by mantle heterogeneity. In Supplementary Fig. 9, we compared our results with the dataset after mantle and ray path corrections were reported most recently 31 where 0 is the Voigt average velocity. As the velocity anisotropy does not show an obvious change with pressure from 330 to 360 GPa, the data calculated at 360 GPa and 6000 K were used for comparison. The calculated travel-time residual of FeH0.25 is able to account for the IC anisotropy with faster velocity along the polar direction. We also calculated velocity anisotropies in the situation of 25%, 50%, and 75% aligned crystals with basal planes parallel to the rotation axis ( Supplementary Fig. 10). In comparison with the recent dataset of seismological observations, approximately 55-70% of LPO is needed to meet the observed residuals of PKPbc-PKIKP, and ~70-80% of LPO is needed to fit the PKPab-PKIKP data. In particular, the IEP model presents milder

Supplementary Discussion 5. H-ion anisotropic diffusion in Fe-H alloys
In this section, we provide detailed methods and discussions for the diffusion properties of Fe-H alloys. The migration barrier energies along different migration paths are calculated using the climbing-image nudged elastic band (CINEB) method 43 Supplementary Fig. 4).
The stable interstitial site is the octahedral site (O), and the tetrahedral site (T) is metastable.
We considered three diffusion paths for H-ions. As shown in Fig. 2a FeH0.25 at 360 GPa and 6000 K, the c/a ratio is 1.626, which is smaller than the ideal value of 1.633.
As the distance of nearby Fe atoms along the c-axis is shorter than the distance in the basal plane, H-ion migration in the basal plane has a stronger interaction with the Fe sublattice. Thus, the c/a ratio is a critical parameter governing the anisotropic diffusion behavior of interstitial impurities in hcp-Fe alloys, and diffusion in the basal plane can be favorable when the c/a ratio is larger than the ideal value. In the case of hcp-FeH0.25, the c/a ratios are lower than the ideal value under IC conditions, and H-ion diffusion along the c-axis is more favorable. From 2000 to 6000 K, the barrier energies for the three paths mostly decrease with temperature, and the barrier energy along the caxis is the lowest (Fig. 2c).
We also used AIMD simulations to study the diffusion anisotropy in hcp-Fe-H alloys. In the AIMD simulations, the diffusion coefficient is defined as: where d is the dimension of the lattice on which ion hopping takes place. Mean square displacement GPa are calculated and plotted as a function of reciprocal temperature (Fig. 2d). The values are fitted with an Arrhenius equation: where ∆ is the activation enthalpy, A is a preexponential factor, k is the Boltzmann constant, and T is the temperature. Based on Eq. 21, ∆ along different paths were calculated, and the activation enthalpy along the c-axis was the lowest (Fig. 4d). This result is consistent with the migration barrier energy calculated using the CINEB method, suggesting that H-ion diffusion along the c-axis is the energetically most favorable path in FeH0. 25. The diffusion coefficient of defects in the lattice can be described by the following equation: where Z is the number of equivalent migration paths, F is the correlation factor, is the migration distance, is the attempt frequency, ∆ and ∆S are the migration enthalpy and entropy, is Boltzmann's constant, and T is the temperature. As shown in Supplementary Fig. 11

Supplementary Discussion 6. H-ion transportation in an external electric field
External electric fields alter ionic transportation behavior in superionic ice 44 Table 3). The total energy in MD simulations was calculated by counting the contributions from potential energy, kinetic energy and electric-field-induced energy change.
The MSDs of the H-ion along the a-axis and c-axis were calculated and are shown in Supplementary   Fig. 14. The MSD changes with the applied electric field, and diffusion along the c-axis is more favorable when the electric field points to the c-axis. However, the electric field along the a-axis does not significantly increase H-ion diffusion along this direction. In particular, when the field intensity is 0.1 V Å -1 , H-ion diffusion along the c-axis is still preferred (Supplementary Fig. 14).